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Abstract 

In this paper we explore the use of an equation of motion decoupling method as an impurity 
solver to be used in conjunction with the dynamical mean field self-consistency condition for the 
solution of lattice models. We benchmark the impurity solver against exact diagonalization, and 



^ ■ apply the method to study the infinite U Hubbard model, the periodic Anderson model and the 



pd model. This simple and numerically efficient approach yields the spectra expected for strongly 
correlated materials, with a quasiparticle peak and a Hubbard band. It works in a large range of 



T3 

parameters, and therefore can be used for the exploration of real materials using LDA+DMFT. 



PACS numbers: 71.27.+a,71.30.+h 



I. INTRODUCTION 

Dynamical mean field theory (DMFT) was developed over the past 15 years into a pow- 
erful tool for the treatment of strongly correlated electron systems^^. DMFT is based on 
the idea of mapping a complicated lattice model onto a single impurity model coupled to 
a noninteracting bath. It relies on the observation that the self energy S(k, iu n ) becomes 
k-independent in infinite dimensions d = oc 4 , making a single site treatment with only 
temporal fluctuation exact in this limit. The DMFT approach derives its strength from the 
fact that it becomes exact in this nontrivial limit of d = oo or infinite lattice coordination. 
Perhaps surprisingly, DMFT proves to be a very good approximation even in d = 3 dimen- 
sions. By replacing complicated models with a single impurity model, the DMFT equations 
can then be solved with one of the methods that have been developed to solve the Anderson 
impurity model. 

The study of correlated materials has until a few years ago been conducted with two ap- 
proaches that are very different in spirit. On the one hand, density functional theory (DFT) 
calculations in the local density approximation (LDA) has proven invaluable in the deter- 
mination of the electronic structure of real materials but there are a number of strongly 
correlated materials where its predictions are even in qualitative disagreement with ex- 
periment. On the other hand, the study of model Hamiltonians has provided qualitative 
understanding of may systems with strong correlations but due to its dependence on pa- 
rameters this method lacks predictive power for new materials. The combination of the two 
approaches in the form of LDA+DMFT 5 promises to deepen our understanding of strongly 
correlated materials as some initial successes demonstrat e - 6 ! 7 ! 8 . 

DMFT reduces the quantum many body problem to a single site problem, namely an 
impurity model in a medium, and requires the solution of an Anderson impurity model 
for arbitrary values of the bath. When the LDA+DMFT is carried out self-consistently 
in a multiband situation, the impurity model has to be solved many times, and becomes 
the bottleneck of the LDA+DMFT algorithm. Therefore it is important to find impurity 
solvers that are reliable and computationally cheap. Currently, the usual choices for solv- 
ing the Anderson impurity model in the framework of LDA+DMFT are quantum Monte 
Carlo (QMC), the non-crossing approximation (NCA), and the iterated perturbation theory 
(IPT). Nevertheless, each of these methods has some drawbacks limiting its range of appli- 



cability. The QMC method is essentially exact, but becomes prohibitively expensive at low 
temperatures and for high interaction strength U. The NCA approximation, applied to the 
impurity model, exceeds the unitarity limit at low temperatures and leads to pathologies in 
the solutions of the DMFT equations. The extension of the IPT scheme, a method which 
was very successful at arbitrary filling in the one orbital situation, has encountered difficul- 
ties in its extension to the multiorbital case. This provides the motivation of this article to 
investigate the usefulness of a previously known decoupling scheme in the context of DMFT. 
The method for the solution of the Anderson impurity model proposed here aims at 
working with an arbitrary noninteracting DOS as input. Nevertheless, we intend to show 
that even for the solution of model Hamiltonians like, e. g. the Hubbard Hamiltonian, a 
DMFT scheme with a closed set of equations gained from a decoupling scheme is superior 
to the direct solution of that Hamiltonian with decoupling methods. 

II. THEORY 

The method of writing equations of motion (EOM) for the Anderson impurity model 
and decoupling them in order to close the system of equations has a long histor y 9 ! 10 ^ 1 . In 
the derivation of the integral equation for the solution of the infinite U Anderson impurity 
model we follow the approach and the decoupling scheme of T. Costi^ 2 -. The Hamiltonian 
for a mixed valent impurity is 12 

H = Y, £ ^inCkn + Y, Sin*™ + E 70*°° + £ (V£ C + X " + V kn X^C kn ) (1) 

kn n kn 

We determine the equation of motion for the m-channel f-electron Green's function (GF) 
F m (u) = ((X 0m ;X m0 )) by writing 

uj((X 0m ; X m0 }) = ([X 0m , X m0 ] + ) + (([X 0m , H];X m0 )) (2) 

and evaluating the appearing (anti-) commutators. The result is 

(w - £ f )((X 0m ; X m0 » =(X 00 + X mm ) + Y Vk(((X°° + X mm )c km ; X m0 » 

k 

(3) 

k 

assuming that the hybridization V k does not depend on the z component m of the angular 
momentum J. The abbreviation ej = Ej m — Ef was introduced. The averages over the X 



operators are (X 00 ) = 1 — n/ and (X mm ) = rif/N where the total number of / electrons is 
calculated as 

n, = ~J&Jf(<J)FW), (4) 

with the notation F m (u) = F^(u) + iF^iu). For the higher order Green's functions on the 
rhs of Eq. Q we also write equations of motion: 



u - e k )(((X w + X mm )c kn , X mU » = (5) 

V k ((X 0m ; X m0 » + J2 V q ((c+ n X 0n c kmi X m0 » + £ V q ((c qn X n0 c km] X m0 » 



q q 



(u-e k )((X nm c kn ,X m0 }}= (6) 

-(X n0 c kn ) + J2 V q ((X n \ m c kn , X m0 )) + J2 V q ((X 0m c+ n c kn ; X m0 » 
i <? 

We now employ the decoupling scheme already given in Ref.— that conserves particle number 
and angular momentum (n ^ m is assumed): 

«c+ X°% m ; X m0 » ~ ( C + X°») ((c km , X m0 » (7) 

((c qn X n0 c km , X m0 }} ~ (c qn X n0 ) ((c km , X m0 }} (8) 

((X n0 c qm c kn ; X m0 » ~ (c fen X"°> «c ?m ; X m0 » (9) 

P^c+c*,; X m0 » ~ ( C + c fen > <(X 0m ; X m0 » (10) 

Note a sign difference in Eq. (jHJ) with respect to Ref.—. The Greens function {{c qm ; X m0 )} 
appearing here can again be determined from its equation of motion 

((c qm , X m0 » = J^((X 0m ; X m0 » (11) 

U) — E q 

This leads to the equation from which the / electron Greens function can be determined: 

F (lu) = l ~ n f + ! W + J iM (u) 

mK ' w-e / -A(w)+J 2 (w)-A(w)Ji(w)' l ; 

with the sums over correlation functions 

J 1 H = -V-^(X n0 c fcn ), (13) 

h{«>) = ~ £ ^^(4nCkn) , (14) 

kg 



and the hybridization function 



AM 



v- n 



z^ 



10 - e k 



(15) 



In the degenerate models we study in this article, the sums over n with n ^ m simply lead 
to factors of N — 1. The sums over k and q can be simplified further. 

To that end, we replace the correlation functions by integrals over the imaginary part of 
the corresponding Greens function 

(4n c gn) = ~~ du'f(u') Im ((c qn ; 4J . (16) 

The conduction electron Greens function {(c qn ; cjj" n )) is determined from its equation of motion 



\C qn ', C kn ) 



>kq 



v k v q 



U- S r , 



(u-e k )(w — e q ) 



((X Un ; X n0 )) 



Now in order to simplify the sums in Eqs. (fTHj) and ()14|). we employ the identity 



{u — e){uj' — s) uj 1 — u [uj — e u' — e 
This allows us to identify occurrences of the hybridization function ()15j) . and we find 



(17) 



(18) 






N - 1 f du>^ 

7T 

N-l 



71 



duj 



uj — uj 

uj' — UJ 



F:( W ')AH-Im{F m ( W ')A( W ')} 



(19) 



-A"(uj') + A(uj)lm{F m {uj')A(uj')} - lm{F m {uj')A(uj') 2 } 



(20) 
Eqs. (|12|l together with (|19jl . (}20|) and the definition of nj (|1]) form an integral equation for 
F m (u)) that can be solved iteratively. In order to compute the integrals of Eqs. (fTHjl and (J20J) 
we introduce the following real functions: 

A m (uj) = -/(w) ImF m (w) 
B(w) = /(w)ImA(w) 

C m H = /Hlm{F m ( W )A( W )} 
D m H = /( W )Im{f m ( W )A( W ) 2 } 



(21) 



Now the integrals read 

7T *-— ' 

7,M = i £ 



A.:,! , ^M + 



n^m 



IT UJ — UJ 

duj' B(uj') 



duj' C m {uj') 



71 UJ — UJ' 



71 UJ — UJ' 



duj' C m (uj') f duj' D m (uj') 

L\(UJ) I 7 + 



71 UJ — UJ 1 



71 UJ — Uj' 



(22) 



Thus, the calculation of the integrals reduces to simple evaluation of Kramers Kronig inte- 
grals. The imaginary part for example of the first such integral is — i7rA m (u'). 

I turns out that this set of equations on the real frequency axis is solved easily for the 
Anderson impurity model, but as we add self consistency conditions in order to solve more 
complicated models in the DMFT approximation, convergence depends strongly on a good 
initial guess of the solution. For this purpose, we write equations analogous to Eqs. (|T2|). 
(|19)L (J20J1 and (@J) on the Matsubara axis. Matsubara Greens functions are much more 
smooth than their counterparts on the real frequency axis and thus converge more easily. 
Nevertheless, the calculation of the Greens function on the imaginary axis does not make 
the real axis calculation redundant: Firstly, the analytic continuation to the real axis is only 
accurate for low frequencies due to a lack of high frequency information in the Matsubara 
Greens function. Secondly, the dependence of the imaginary frequency grid on temperature 

iu n = (2n + 1)ttT (23) 

means that at high temperatures, the low frequency part of the Greens function is very badly 
resolved, while at very low temperatures, an inordinate number of imaginary frequencies is 
necessary to describe the Greens function for all frequencies for which it significantly differs 
from zero. This means that from a practical point of view, the Matsubara Greens function 
is best calculated at an intermediate temperature, providing via analytic continuation a 
sufficiently accurate initial guess for the iterative solution of Eq. 112(1 on the real axis. 
This problem of the Matsubara formulation is not related to the well known difficulty in 
performing analytic continuation to the real axis. 

All equations of motion are almost unchanged when we go over to Matsubara frequency 
iu n , e. g. Eq. (|T7|) becomes 

<fe„; <,»,„. = j-%- + —W±-—AX°«; .Y»»»,„„ . (24) 



mo,, — e 



i 



(iu n -e k )(iuj n -e q ) 



Correlations have to be calculated as 

(ctnCkn) =TJ2 «<*»; C+ ))^e^° + (25) 

iul n 

which replaces Eq. (jlfij) for that purpose. In order to simplify the equations, we employ the 

analog of Eq. (fT£J) , namely 

V 2 1 f V 2 V 2 } 

Vk J v h v k i (26) 



(iuj n - e k ) (iu' n - e k ) iuj' n - iuj n { iu n - e k iu' n - e k 



and we identify occurrences of the hybridization function 



A(iu n ) = V , V " . (27) 

*-£ ioj n - e k 



This leads to the system of equations 



F m {iu n ) = l-n f +^ + S^) 

iu n - Sf - A(iu n )(l + Si[iu n )) + S 2 {ico n ) 



with 



&(*"») = T E A( T; ) :f^" ) i 1 + A(-;)^«)}e-" 0+ (30) 

"I. 

n f = NTY,F m (^'n)e l<0+ . (31) 

With the replacement 

T^A^)e^ 0+ — ► -- f doo'f(uj')lmK(uj') (32) 

we can easily recover Eqs. (fTT?|) and (|20j) from (|29jl an d (j3*Uj) . It is important to note that good 
convergence of the self-consistent solution of the system of equations depends crucially on the 
proper treatment of the slowly decaying high frequency tails of the addends of Eqs. (|29 p -([31| ) . 
A high frequency expansion of these addends was performed to determine the coefficients of 
the terms proportional to l/iu> n and l/{iuj n ) 2 . These terms were subtracted from the sums, 
and their value was determined analytically. 

A. Hubbard model 

We now proceed to investigate the usefulness of the impurity solver detailed above in its 
application in the DMFT context. Our application of the method to three lattice models 
is an exploratory study concentrating on a small number of important properties only. It 
is not the intention of this article to go into detail for each of the three models. We first 
investigate the Hubbard model in order to study the quasiparticle scaling of the Hubbard 
band with degeneracy N . 



We consider the Hubbard Hamiltonian 



ija ia iacr' 

a^a' 

where the spin and orbital index a runs from 1 to N. For this model, we have to solve the 
AIM with the self-consistency condition 

A(iu n ) = t 2 G ooa {iu n ) . (34) 

For the derivation of this equation, see App. |XJ 

B. Anderson lattice 

We study the application of the U — oo impurity solver to the Anderson lattice in order to 
learn how this new approach compares to the straightforward decoupling of the equations of 
motion for the periodic Anderson model 12 . We consider the periodic Anderson Hamiltonian 

#pam = -J2 tc ij c t c i° - Yl tij&fj* +£ cYl c ^ Ci ° +e fYl f*f*° 

ija ija ia ia 

' / J ^yia^iaJia + *iaJia C iv) ' 7f / 4 J ioJ iaj i a i 1 ia' 

ia iaa' 

a^a' 

In this case, the self consistency condition for the / electron Greens function is 

Gf*(iu n ) = fds p (e) (iu n + fi-e° f - E f (iu n ) - — ^ ) (36) 

Here, the self energy is determined from the equation 

Go 1 fan) = G^ x {ioj n ) + E(iu n ) (37) 

and the Weiss function Q Q (iu n ) is related to the hybridization function A{iuj n ) by 

Go 1 fan) =iu n + fi-e° f - A(iu n ) (38) 

The derivation of these equations is contained in App. |Bj 



C. pd model 



In order to study the Mott transition with the U = oo impurity solver described above, 
we consider the Hamiltonian 14 

H = ~ J2 V v [ d tPj<r + v]A° + £ pJ2 p t- p i° + £d Yl d ^ di ° + Ud J2 d t] d ^ d ti d H • (39) 

ijcr jtr ia i 

This Hamiltonian, which we call pd model here, is similar to the Anderson lattice Hamil- 
tonian if the conduction electron dispersion is taken to be a constant e^ = e p and if the k 
dependence of the hybridization Vk is retained. This changes the local conduction electron 
Greens function: 



= / de pp d (e) (iu n + p-e d - S rf (iw 
= (p / dep pd (s)(( p ( d -e 2 j 
G l ° caX (ioj n ) = / dep pd (e) (iu n + p-e p - - — 
= Cd I dep pd (e)(CpCd-£ 2 J 



iu n + p - e p 



+ p - e d - Y, d (iu ri 



-i 



(40) 



with the abbreviations ( p = iui n + p — e p and ( d = iu n + p — e d — Ti d (iui n ). Here, p pd (e) 
stands for the density of states associated with the hybridization V&. Noting that 



de 



Ppd{e) = l_ 
x 2 — e 2 2x 



de 



P P d{e) 

x — e 



ds 



P P d{e) 



D{x) 



x 



x + e 
for symmetric p pd (e) 



j-(D(x) - D(-x)) 



(41) 



we can write G l ° cal (iuj n ) and G^ ocal (iu; n ) as Hilbert transforms 



G i r\iu n ) 

We use a semicircular form for p pd (e): 

Ppd(e) = 



Kdf, 
C P 






(42) 



2*&V 4 ft--" J 



(43) 



where t pd is the strength of the hybridization between p and d levels. 



It is worth pointing out, that this method reproduces an important aspect of the exact 
solution of the DMFT equations within the context of the pd model. Namely, it produces 
a first order phase transition between a metallic and an insulating phase, which is manifest 
by the existence of two DMFT solutions for the same range of parameters. 

III. RESULTS 

A. Hubbard model 

First of all we test the performance of our impurity solver by comparing it with the results 
of exact diagonalization (ED). For this purpose, we employ the code published accompanying 
the review of the DMFT method in Ref.— , modified to U = oo. The Hubbard model is 
solved in the DMFT approximation. The self-consistency condition for the Hubbard model 
is realized by minimizing the function /(£&, V&) = J2 n \t 2 G(iu n ) — J2k Vic/i^n ~ £ k)\ with 
respect to the parameters Ek and Vk ■ Here, the exact diagonalization has been performed 
with N s = 6 sites which are divided into 1 site for the impurity and 5 sites for the bath. 
Thus, the hybridization function A(iu n ) is represented with 5 poles. This leads to a finite 
number of poles instead of a smooth function in the spectral function as well. Fig. ^ shows 
the comparison of the densities of f electrons as a function of the impurity position Ef 
(which is related to the chemical potential by \x = —ef)- The comparison shows that at 
high temperature T = 0.5, the results of exact diagonalization and EOM are virtually 
indistinguishable while for a lower temperature T = 0.03, the densities differ slightly for 
impurity positions between -1 and 1. 

In Fig. EJ we compare the imaginary parts of the Greens function for a density of nf = 
0.84. The slight differences in the n/ versus \x curves of Fig. Q mean that this density 
is achieved for p = 0.6 in the case of ED and for p = 0.53 in the case of EOM. The 
imaginary parts of the Greens function on the Matsubara axis shown in the inset are very 
similar. Thus, the main figure shows the more demanding comparison of the densities of 
state. The continuous line represents the DOS from the EOM method gained by analytic 
continuation in the Pade approximation, while the long dashes stand for the EOM result 
on the real axis. The dashed curve with the five poles is the result of DMFT on the basis 
of exact diagonalization. The figure shows that the distribution of spectral weight between 



10 



the Kondo peak around the Fermi level at u = and the Hubbard band is similar in both 
methods, but the EOM method leads to a better overall shape of the spectral function. We 
conclude that the EOM method results compare well with ED, giving us confidence that it is 
a useful approximation. Even for this low number of N s = 6 sites, the exact diagonalization 
requires an order of magnitude more CPU time than the EOM method. 

Fig. El shows the carrier density as a function of the impurity position £f. The impurity 
position corresponds to the chemical potential, only with opposite sign /x = —£f. Due to 
the infinite interaction, the maximum filling is 1 electron per site. In other words, the 
upper Hubbard band that could hold a second electron at finite U has been pushed to 
infinite energy. While at low temperature T = 0.03 the nj versus /i curves at different 
degeneracies N = 2 to N = 14 nearly coincide (see Fig. 131(b)), they differ considerably at 
high temperature T = 0.5 (see Fig. El (a)). 

In Fig. |U we show examples of the spectral function for degeneracies between N = 2 and 
A = 14 for high and low temperature. While at T = 0.5 the spectral function is nearly 
unstructured, at T = 0.003 a broad Hubbard band and a quasiparticle resonance at zero 
frequency u = can be distinguished. The weight of the Hubbard band diminishes as 1/N as 
the degeneracy A increases while the intensity of the Kondo peaks remains nearly constant. 
Note that the spectral functions in Fig. |U resulting from the DMFT self-consistency contain 
no spurious side bands as those calculated by directly decoupling the equations of motion 
produced by the Hubbard Hamiltoniani 3 . In our calculation, the imaginary part of the 
Greens function outside the Hubbard band and resonance is exactly zero. 

B. Anderson lattice 

Fig. 03 shows examples for the conduction electron and the strongly correlated / electron 
spectral functions (dashed and full lines, respectively). In Fig. (a), the hybridization 
between the two bands is small (V 2 = 0.01) while in Fig. El (b) it is rather large (V 2 = 0.2). 
Correspondingly, the conduction electron DOS shows only a small dip at the position of the 
f band for a low value of the hybridization. Interestingly, we find a Kondo resonance at the 
Fermi level in the / electron DOS. This resonance was absent in the decoupling approach 
to the periodic Anderson model of Ref.— . 
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C. pd model 

We investigate the pd model Hamiltonian Eq. (|39j) as a function of the separation Ao = 
e p — Ed and of the hybridization strength t p d between the two bands. From the analysis in 
Rcf. 11 of the finite U version of this model, wc expect a metal insulator transition to occur 
at a fixed density n to t = 1 if we vary the level separation Ao at a given t p d- Fig. El (a) a 
shows the result of this calculation at a fixed t p d — 1. The temperature was taken to be 
T = 0.01. For level separations Ao = and Ao = 0.5 , the density n tot around n tot = 1 
changes smoothly as a function of chemical potential. But beginning with A = 1.0, a 
charge transfer gap g\ = fJ,(n to t = 1 + ) — Mn toi = I - ) begins to open up. Thus, the physics 
discussed in Ref.— for finite values of U can be also found for U — oo. The critical value 
at tpd = 1 is Ao = 1. Note that the Ao = 4 U = 8 result in Fig. 1 of Ref.— compares well 
with the A = 4 U = oo curve of this work's Fig. El (a). If we increase the hybridization 
strength to t p d = 4 (see Fig. El (b)), we find that the critical A for the metal to charge 
transfer insulator increases to Ao ~ 4. In Fig. El (a), we also note the transition at a total 
density n tot = 2 from a metal at higher level separation A to a band insulator with a gap 
g 2 = n{n tot = 2 + ) — nijitot = 2~) . For the higher value of the hybridization strength t p d , 
the system is a band insulator at n tot = 2 for all studied level separations Ao- 

An important question in the metal to insulator transition of Fig. El concerns the existence 
of a coexistence region. We can show that such a coexistence is indeed found with our 
method. Fig. shows spectral functions for d and p electrons at a hybridization strength 
tpd = 1, a separation A = £ p — Ed = 1 and a chemical potential fi — Ed = 0.3. The calculation 
was performed for a temperature of T = 10~ 5 , and care was taken to resolve the sharp peak 
of the noninteracting Greens function G p (u) at e p = 0.5 with the help of a logarithmic 
grid. The full line shows the converged result of a direct calculation at \x — Ed = 0.3. A 
quasiparticle peak at u = for both the correlated and the uncorrelated electrons makes 
this a metallic solution. The dashed line was obtained by starting the calculation in the 
insulating region at /i — Ed = 0.5 and lowering the chemical potential in steps of 0.01. At 
/i — Ed = 0.3, the solution is still insulating as no quasiparticle peak has formed. 
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IV. SUMMARY 

A method to solve the Anderson impurity model with the help of equations of motion 
and decoupling has been tested for its suitability as an impurity solver in the framework of 
dynamical mean field theory. The application to three lattice models in infinite dimensions 
and for infinite interaction strength U shows very encouraging results. In the application to 
the Hubbard model, we see a correct quasiparticle scaling of the Hubbard band with respect 
to the degeneracy. In the periodic Anderson model, we find a Kondo resonance which is 
absent in a direct decoupling of the equations of motion. This underlines the usefulness 
of the approach chosen here: To use a decoupling scheme for the solution of the Anderson 
impurity model which is then employed to solve lattice models in the DFMT approximation. 
Interestingly, the application of our approach to the pd model yields a coexistence of the 
insulating and metallic phases. The extension of the U — oo approach discussed here to 
finite values of the interaction strength U is possible and in preparation. The numerical 
efficiency of the method makes an application in an LDA+DMFT context feasible. 
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APPENDIX A: DMFT SELF CONSISTENCY CONDITION FOR THE HUB- 
BARD MODEL 

The partition function corresponding to the Hamiltonian of Eq. (J2SJ) is 

Z= I ' Y[Vc ia Vc la e- s (Al) 

** i 

with the action 

rP q r0 r 

u icr </ L ija j CT 

(A2) 



+ -z- ^ Cjff(r)ciff(r)ciff'(r)ciff'(r 
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where the Fermion operators c[ a , q ct of the Hamiltonian have been replaced by Grassmann 
variables Ci a (r), Ci a {r). 

The cavity method now requires that we focus on one site % = o and separate the Hamil- 
tonian (J33j) into three parts, one relating to site o only, one connecting this site to the lattice 
and one for the lattice with site o removed: 



H = H + H C + # (o) 



H n 



H r 



U 



/ , C oa C oa + „ £ , 



C ocr C ocrC oa iC oa i 



n 



CFU 



^io^icr^oa ~T~ <■ oi^ OIJ C-i<j 



H (o) 



u 



/ , ^ijC^Cju fi y j c ia Ci a + y j c ia Ci a c ia , 



Cirr' 



(A3) 
(A4) 

(A5) 
(A6) 



i^OJj^OCT 



i^ou 






The three parts of the Hamiltonian correspond to the action S of site o, the action AS for 
the interaction between site o and the lattice, and the action S^ of the lattice without site 



S = dr ^2c oa {r)^- // 

AS = - dr 2]t i0 Ci ff ( 



CoAt) 



-E 

crj^cr' 



Coa(T C oa (T C oa >(T C oa >(T 



■i=£oa 

u 



T)C oa {T) + t oi C oa (T)c ia (T) 



i]Cia\' )Cjcr\' ) 



(A7) 

(A8) 
(A9) 



+ — 22 c ^( r ) c ^( r ) c ^'( r ) c ^'( r ) 



The aim is now to integrate out all lattice degrees of freedom except those of site o in 
order to find the effective dynamics at site o. In that process, the action S remains un- 
changed, the terms of AS are expanded in terms of the hopping t which becomes small 
with increasing dimension and averaged with respect to the action S^°\ Defining AS(r) via 
AS = J Q dr AS(t) the partition function is 



Z=f Vc OCT Vc oa e- s ° J J] Vc m Vc la e- s{o) e~ H ** S V 



(A10) 
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Now we can expand the last exponential function as 

e -J drAS(r) = 1- drAS(r) + - / dry \ dr 2 AS(n)AS( 
Jo 2! J J 



(All) 



(A13) 



Taking into account that in general an operator average with respect to an action S can be 
expressed as 

{A)g = JU^V^^,^] = z -iJY[Vc a Vc a e- s A[c a ,c a ] (A12) 

we can consider the second functional integral in (JA10J) to average the terms of the expansion 
(|A11|) with respect to the lattice action S^: 

Z= I Y[Vc 0ff Vc 0(T e- s °Z s(o) h - J dT(AS(r)) s(o) 

+ 1 f dn f dT 2 (AS(n)AS{T 2 )) sW 

^- Jo Jo 

Here, the partition function of the lattice without site o is abbreviated as 

Z 8 M = / II v c a Vc a e- s(o) . (A14) 

•* i 

Now the terms in (JA13J) with odd powers of AS will average to zero. For example, 

(AS(t)) S (o) = ^t io (c iCT (r)) 5 ( o) c OCT (r) +t oi c OCT (r)(c iff (r)) 5(o ) = 0, (A15) 

ia 

because the average ()s(.°) ac ts on all sites except o. The next average in (|A13j) yields 
(AS'(ri)AS'(r2)) s .(o) = (T T ^^ Cj (T (ri)c 0(T (r 1 ) +t oi c 0<T (T 1 )c i(T (ri) 

ia 

^t jo c jV (r 2 )c OCT /(r 2 ) + t jC 0ff /(r 2 )c JCT /(r 2 ) 

ja' 

^2 ^ojC 0ff (r 1 )(T T c iCT (ri)c JCT /(r 2 )) 5(o ) c OCT /(r 2 ) + ^ ^ioCo<7(Ti)(r T c iCT (r 1 )c jff /(r 2 )) s .( ) c OCT /(r 2 ) 

tj(TCT f ijacr' 

2 ^2 t ^ jC OCT (r 1 )(T T c iCT (r 1 )c J v(r 2 )) s( o) c OCT /(r 2 ) 

ijacr' 

2 } J t io t oj C off {Ti) {T T c iff {Ti)c jff {T 2 )) S (°) c ou {t 2 ) 

ija 
-2 ^ tiotojCo<r(Tl)G < l°> (T (Ti - T 2 )c 0(7 (t 2 ) 

(A16) 



S(o) 



IJV 
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The imaginary time ordering operator T T enters because the path integral leads to imaginary 
time ordering. Only terms with a = a' contribute as we are considering a paramagnetic state 
and thus {T T Ci a {Ti)Cj a i{j 2 )) S ( 0) = 5 - -/(T T Cj -(r 1 )cj -(r 2 ))5(o). We have identified the average 
with the cavity Greens function G^ ff {r\ — r 2 ) = — (T r Cj cr (r 1 )c+ 7 (r2)) 5 (o), i. e. the Greens 
function of the Hubbard model without the site o. Now we have for the partition function 

5)X 



Z= \{Vc OCT Vc OCT e- s °Z s{ o 

f f 13 f 13 (o) 1 ( A17 ) 

x <M - / dn dr 2 ^ tj t jC 0CT (ri)c 0CT (r 2 )G ! -° ) CT (ri - r 2 ) + . . . > 

We would like to write the bracket {} in (|A17|) again as an exponential function in order to 
identify an effective action S e s. 

Z = J II Vc oa Vc oa e- s * a (A18) 



Noting that the next term in the expansion of (jA17|) would read 
,/3 r p r p r p 

. _ i x 

u *2 ii h o 

x tiioti 2 otoj 1 t j 2 G ili2 j 1 j 2(T (Ti r 3 , r 2 r 4 J 

we can write for the partition function (JA17J) 



nP ,0 r-P rf 

I dn / dr 2 / dr 3 / dr A 2J c oct (ti)c oct (t 3 )c 0(T (t 2 )c oct (t4) 
>>o Jo Jo Jo 



(A19) 



»)X 

rP r P 



Z= l\\Vc oa Vc oa e~ s °Z s{o) 
J i 

r °° rP rP 

x exp<^ - ^ Z^2 / rfr i • • • / rfr 2n c OCT (ri) . . . c OCT (r 2n _ 1 )c OCT (r 2 ) . . . c OCT (r 2 „) x , A2Q v 

/ , *no-- -tinotoh ■ ■ ■ t oj n G iiAnji<jnCr (T 1 . . .r 2n _i,r 2 . . .r 2n ) > 



jl,—,jn 



All terms but the first in this sum over n turn out to be at least of order 1/d so that they 

vanish in the limit of infinite dimension d — 00. Thus, in this limit we find for the effective 

action 

rP rP 

S c s = So + 22 / dn I dr 2 c oct (ti)c oct (t 2 ) V, t io t oj Glf a (ri - r 2 ) 



/ dr y^c OCT (r)(- A i )c OCT (r) + — V" c 0ct (t)c 0(T (t)c 0(T /(t)c 0ct '(7) 

cryV 

/•/3 rP 

+ 22 dT l dT2C OCT (T 1 )c oa (T 2 )y2t io t oj G\°l(T 1 -T 2 ) 

„ Jo Jo 



(A21) 
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van 



tA&nn -J2 tiot °i G i%( Tl ~ T ^> 



(A22) 



u 



and introducing the Weiss field 

we finally get 

■S'eff = ~ /J / dri / ^r 2 c 0ff (ri)^ 1 (ri-r 2 )c OCT (r 2 )+ / rfr — ^ c 0(T (t)c 0(T (t)c 0(T /(t)cw(t) 
„ Jo Jo Jo l , 



u+a' 



The equation 



ij a *~*ijcr ( - Jr io a^ oo a oj & 



(A23) 



(A24) 



is needed to relate the cavity Greens function to the Greens function of the lattice G{j a . 
Going from imaginary time to imaginary frequency and combining with (|A24|) . the Weiss 
function (jA22|) reads 



g a 1 (iu n ) =iuj n + fi-^2 tiotojGif ff (iu n ) 

•>n + A 4 — / J tjotpj Gij a (iuj n ) — Gi 0a (iu) n )G ooa (iuj n )G j a (iuj n ) 



lUJ r , 



(A25) 



If we now go from real space to k space we can simplify this equation. Introducing the 
Fourier transform Gka via 



G ij(T (iu n ) =^e ikR ^G ka {iu n ) 



(A26) 



we find 



y^JioGioajiUn) = y^Ji o y^ j e lkRio G k(T (iuj n ) = ^e k G ka {iu n ) 

i i k k 

/ jtjotojGjj cr(iuj n ) = y tj t j 2 J e - 1 ^Gkuii^n) (A27) 

ij ij k 

k i j k 

In the general form of the Greens function G^ a (iuj n ) = iu n + fj, — e k — S cr (zu; n ) we introduce 
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the abbreviation £ = iuj n + ji — "E ff (iu n ) to get G k l(iuj n ) = £ — e k and determine the sums 



22e k G k(T (iu; n ) = \] 



^A: 



^ - e + e 



•i+E 



"1 + f E G k<r( iuJ n) = -1 + £,G ooa {iu n ) 



£ 



£-£fc 



J~]elG ka (iv ri 



it s fe 



^fc(^fc ~ + £ fc£ 



(A28) 



E^E 



£fc 



A- 



^-^fc 



- £(-1 + £G 00(T (^n)) = -£ + eGooA 
With this, the Weiss function (JA25J) becomes 

Q~ l {iu n ) = iu n + [i - s ^ j e\G k(y {iu n ) + ( ^2e k G k(T (iu n )) G~^(iu n ) 

k k 

= iu n + [i + i - £ 2 G 00(J {iu) n ) + (-1 + iG ooa {iu) n )) (-G~^{iu n ) + f) 
= zw n + // - £ + G-l a (iu n ) = E (iu n ) + G^ a (iu n ) 



(A29) 



This equation G^ a (iuj n ) = Q~ 1 {iui r ^) — E ff (wi> R ) is the Dyson equation for the local Greens 
function. 

The effective action (|A23|) can now be interpreted in terms of the Anderson impurity 
model, i. e. the Anderson impurity model gives rise to an action which becomes identical 
to (jA23|) if an additional self consistency condition is fulfilled. The Hamiltonian for the 
Anderson impurity model is 

h = y1 *#*><** + E &*&& + Vkfi***) - E vtif' + % E ttfotiu ( A3 °) 



ka 



era 



where a runs from 1 to the degeneracy N. The action corresponding to this Hamiltonian 
will consist of a purely local part S concerning only the / electrons 



S ° = [ dr ^2Ur){^-^)Ur) + ^Y, Ur)Ur)Ur)fAr 






(A31) 



and a part involving conduction band electrons that can be integrated out: 



r0 r . 

S = S + ^E ^( r )( 
Jo k A v 



d 
Cka(r) ( -K- + e k )c k<T (r) + V k c ka (r)f a (r) + V k f ff {r)c k<T {r) 



(A32) 
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Now the partition function for the Hamiltonian (|A30|) is 

Z = J Vf.Vf, J TJ Vc la Vc la e- s = j Vj.Vf, e" 5 " / TJ Vc„Vc 

'*■£ 

ka 

Vf a Vf a e~ s °l[det(— + e k 



x exp< 



d 

c k a(r) ( — + e k )c fe0 .(r) + V k c ka {r) f a (r) + V k * / CT (r)c fc(J (r) 



(A33) 



x 



rf 3 ri 3 / Q \ -i 



u 



u 



&n 



In the last step, the terms involving / electrons V k f a (r) and V k f a {r) were taken as source 
terms, which makes the term in the exponent a Gaussian integral that can be evaluated 
directly. The determinant constitutes a constant factor in the partition function that doesn't 
concern us here. We are left with an action for the / electrons that reads 



rP rP _ / d 

Sf= l dTl l ^ 2 £ / * (Tl) fe-' i 

+ / d^Y.Ur)Ur)JAr)fAr 

■J / 



5> 



,1 i^ + ek) *r 



Ur 2 ) 



(A34) 



(JCF 



If we now compare this to the effective action of the Hubbard model f!A23|) . we see that they 
are identical if we require that the Weiss function Q{j\ — r 2 ) fulfils the condition 

d n-1 



■(- 



G 1 (ti-t 2 ) = -[-£■ -P) S nn + ^2\V k ^[^-+ £k) 8 nr ., 



(A35) 



Going from imaginary time to imaginary frequency, this equation reads 



Q 1 {iu n ) = iu n + n - y^ 



IT/ I 2 



iu n — e k 



(A36) 



Here we can identify the usual definition of the hybridization function A(iu n ) in the Ander- 
son impurity model 



A(iu n ) = ^ 



IT/ I 2 
\Vk\ 



iu n - e k 



(A37) 

If we now equate Weiss functions (|A29|) and (|A36|) we find the DMFT self-consistency 
condition in terms of a prescription for A(iuj n ) 



A(iu n ) = iu n + /j, - T; ff (iu n ) - G ooa (i(jj n ) 



(A38) 
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On the Bethe lattice and with a half band width of D = 2£, we have a noninteracting density 

of states 

1 



Po(s) = ^V4t^^ (A39) 



and thus we can write for the local Greens function 



G ooa {iu n ) = E G k{^n) = E 7— — with f = iu n + y, - Y, a (iu n ) 

k k / , (A40) 

From this we gain the expression 

t 2 G ooa (icu n ) - £ + G- \(tu n ) = , (A41) 

which combined with Eq. (JA38)) leads to a simplified form of the self-consistency condition 

A{iuj n ) = t 2 G ooa {iu n ) . (A42) 

APPENDIX B: DMFT SELF CONSISTENCY CONDITION FOR THE ANDER- 
SON LATTICE 

We again focus on one site i = o and split the Hamiltonian into three parts: 

H PAM =H + H C + H^ (Bl) 

-"o ^c / j C oa C oa ~\- Ef y J oa Joa 

(7 (7 

+ J2( v °° c oJo« + v:jtc oa ) + ^ e foJo«fi>fo*> (B2) 

a ucr' 

cr^cr' 

H c = -H [ tC ie C trCo„ + t C m ct a C ta \ (B3) 

ia 

H {o) = - E ttjC+c^ + e^ctcia + ef^ftM 

+ E ( v ^ c tu + k&c*) + 1 E fcftofof* ( B4 ) 

H c has the same form as in the Hubbard model, but the local part H is more complicated 
as it contains two species of electrons, conduction and / electrons. Nevertheless, we can 
proceed completely along the lines detailed for the Hubbard model above, expanding the 
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action AS arising from H c in order to arrive at an effective action for site o. In this case we 
have 



So =J dr Y,Ur)^ + e f )Ur) + ^Y,fo,{r)f oa {r)f oal {r)f oa ,{T) 

it air' 

+ J2 Co*(t) (tt- + £c)c 0(T (t) + J2( V o^oA^)foa(r) + V *J 0(T (r)c 0(7 (r)) 



and 



,/3 r (3 



f f 

S c s = S + dn / dr 2 c oa (T 1 )c oa (T 2 )y2t io t oj Gly a (T 1 - r 2 

Jo Jo ija 



In the d — > oo limit, the Green's function becomes 



G 1 (iu n ,k) 



ioo n + fi — e°r — T,f {ioj r , 



V, 



v k 

iu n + (x- e k 



Inverting the matrix according to 

'a B 



we find 



M 



G(iu n , k) 



1 



, , C -B 

rxM- 1 

B C I detM I -B A 



iu n + 11 - e° f - Y,f(iuj n ))(iuj n + 11 - e k ) - V k 2 

iu n + Li- e k —V k 

-V k iuj n + fi-e° f - T,f(iu r , 



x 



Thus, we find for the / electron Green's function 

Gf(iu„, k) = liu n + pL - e° f - T,f(iu n ) 
and for the conduction band Green's function 
G c (iuj n , k) = I iu n + 11 - e k 



Vi 



iuj n + 11 - e k 



V k 2 



-1 



(B5) 



(B6) 



(B7) 



fB8) 



(B9) 



(BIO) 



(Bll) 



iu n + /i - e° f - T,f(iu; r , 
We get the local propagators as Gf(R = 0, iuj n ) = ^ k Gf(iu n , k)e lfe ^ =0 ^ by summation over 



k: 



G l p l (tu n ) = Y,G f (tu Jn ,k) 

de p (e) (iu n + /j,-e° f - T, f (iu> n ) 
G l ° cal (icv n ) = / de p (e) (icv n + p-e- 



V(e) 



iu n + ii — e / 

2 \ -1 



(B12) 



iuj n + fi — e°r — Tif(iu}. 
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For computational purposes it is useful to note that for the case of an energy independent 
V(e) = V, G l ° ca " l (iiu n ) can be written as a Hilbert transform D(Q = f^dej^ 



-OO £— £ 

V 2 



G l r\iun) =D(tu n + fi- —— ' ) (B13) 

Rewriting G l j cal (iu) n ), we can likewise reduce the energy integral to the calculation of a 
Hilbert transform: 



G l p\iuj n ) = j ' de p (e)| 



1 



+ 



iuj n + ji — £°r — T,f(iuj n ) 

V 2 1 



,2 ■ 



(B14) 



(iuj n + fi-e° f - S/(iw«)) iw n + n-e- iWn+ ^ e y Sf(ilA)n) 
and with Eq. (jEH 

G " al( ^ } = iu + u-J>-Z f (iu;) + -T-L r y r ^ G l r\^n) (B15) 

If we now assume a semicircular DOS D{e) = 2^2 v^t 2 — £ 2 for the hybridization Vk we can 
explicitly write for the Hilbert transform 

D(Q = ^f2J J £ \_ £ = ^(C-sgn(ReC)VC^4^) (B16) 

Thus, on the Bethe lattice the self consistency condition can be calculated without an integral 
over energies. We also need the Dyson equation 

g Q \iu n ) = Gj\iu n ) + £(iw„) (B17) 

From the high frequency limit of this equation we can find the form of the Weiss function 
0Q 1 (iuj n ) by comparing the terms of the expansion order by order. Expanding Eq. (JB12J) we 

find 

1 / 1 \ 2 

Gf(iu n ) ~ h (e °f — [I + T,f(iu n )) ( ) for iu n — > oo (B18) 

iu n ' \%ijj n ) 

Expanding the inverse, we find 

Gj (iu> n ) « iuj n + /i — e°j- — Tjf{iuj n ) for ito n — > oo (B19) 

Thus, we find from Eq. (|B17|) the high frequency form of the Weiss function: 

g^iiu^^iUn + ^-e) (B20) 
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The hybridization function A(iu n ) contains what we have neglected in the high frequency 
expansion: 

g G l {iuj n ) = iu n + /i - e° f - A(iu n ) (B21) 
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FIG. 1: Density of / electrons as a function of the chemical potential /x = — £/ for exact diagonal- 
ization and the equation of motion method in comparison. The energy unit is the half band width 
D. For the higher temperature T = 0.5, the two methods agree extremely well, while for the lower 
temperature T = 0.03, the exact diagonalization gives slightly lower densities at the same chemical 
potential fi = — £/. Exact diagonalization was performed with 6 sites, and the Hubbard model was 
solved in the DMFT approximation. 
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FIG. 2: Spectral functions (main figure) and imaginary part of the Greens function (inset) from 
exact diagonalization and the equation of motion method in comparison. The temperature is 
T = 0.03, the density of / electrons is rif = 0.84 for both methods. The two methods compare 
well, considering that the exact diagonalization with 5 bath sites has only limited resolution on the 
real axis. 
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FIG. 3: Density nj of / electrons as a function of the chemical potenital /i for the infinite U 
Hubbard model. Energy is measured in units of half band width D. (a) At high temperature 
T = 0.5, nf(sf) differs for different values of the degeneracy N. (b) At low temperature T = 0.03, 
the nf(ef) for different iV nearly coincide. 
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FIG. 4: Density of states Pf(u>) of / electrons for the infinite U Hubbard model. The energy unit 
is the half band width D. (a) At high temperature T = 0.5, there is no quasiparticle resonance at 
ui = 0. (b) At low temperature T = 0.03, the quasiparticle resonance at u> = is clearly developed. 
The weight of the Hubbard band is proportional to 1/N. 
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FIG. 5: Densities of states Pf{oo) and p c {u>) of / and conduction electrons for the infinite U 
periodic Anderson model, (a) At a low hybridization V 2 = 0.01, the / electron Greens function 
is maily a peak at the impurity position; there is no quasiparticle resonance at u = 0. (b) At a 
high hybridization V 2 = 0.2, the hubbard band of the / electron Greens function is split into two 
peaks, and the quasiparticle resonance at uj = is clearly developed. 
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FIG. 6: Total density of electrons as a function of chemical potential for the pd model Hamiltonian 
((HHl)- The hybridization strengths were (a) t p d = 1 and (b) t p d = 4. The plateaus at densities 
ntotai = 1 arid ntotai = 2 correspond to the insulating phases. 
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FIG. 7: Spectral function for the pd model Hamiltonian (|39|) showing the coexistence of a metallic 
and an insulating phase. The hybridization strength was t p d = 1, the pd separation Ao = 1. (a) 
shows the correlated d spectral function, (b) the uncorrelated p spectral function. The full line 
corresponds to a metallic solution, the dashed line to an insulating solution. In (b) the sharp peak 
of the noninteracting DOS at e p = 0.5 is not shown. 
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